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We study the current and shot noise in a linear array of metallic nanoparticles taking explicitly 
into consideration their discrete electronic spectra. Phonon assisted tunneling and dissipative effects 
on single nanoparticles are incorporated as well. The capacitance matrix which determines the 
classical Coulomb interaction within the capacitance model is calculated numerically from a realistic 
geometry. A Monte Carlo algorithm which self-adapts to the size of the system allows us to simulate 
the single-electron transport properties within a semiclassical framework. We present several effects 
that are related to the geometry and the one-electron level spacing like e.g. a negative differential 
conductance (NDC) effect. Consequently these effects are designable by the choice of the size and 
arrangement of the nanoparticles. 

PACS numbers: 73.63.Kv, 73.23.Hk, 85.35.Gv, 72.70.+m 



Programmable self-assembly is one of the most promis- 
ing bottom-up approaches to the synthesis of nano- 
electronic devices. For this technique a template is 
needed which determines the desired device shape. It 
has turned out that DNA is ideal for this purposed be- 
cause of the highly selective binding of single-stranded 
DNA to another strand with complementary bases, it 
can be conformed to a variety of geometrical shapes^. 
By decorating the DNA with metal nanoparticles con- 
ducting material can be created. The attachment of gold 
nanoparticles with a selective size to surfaces^, certain 
biopolymers^ and to DNA^ has already been achieved. 
In this paper we focus on linear arrays of nanoparticles - 
also called (quantum) dots. This geometry is interesting 
for applications e.g. since its electron transport proper- 
ties are robust against unintentional fluctuating back- 
ground charges^. On the other hand, fundamental trans- 
port phenomena can be observed in such systems^. In 
our case, the nanoparticles are attached to a DNA strand 
via a functionalizing ligand shell, see Fig. ^ f° r a wire- 
frame. The DNA strand, which resides on a dielectric 
substrate, is stretched between two macroscopic metal 
leads, serving as electron reservoirs, to which different 
potentials can be applied. A third metal lead serves as a 
gate. 

We determine the current originating from the tunnel- 
ing of single electrons and the corresponding shot noise 
as functions of various parameters like e.g. the applied 
gate and bias voltage, the temperature and the strength 
of dissipative effects. Especially the geometry of the ar- 
ray plays a major role since it is intended to design the 
electron transport properties by controlling the shape of 
the device. From existing studies of arrays of metallic 
nanoparticles with small capacitances and high junction 
resistanceai n *ii*i2ii&i£ii£ii& it is known that such systems 
exhibit nonlinear current-voltage characteristics - also 
called IV characteristics - which can be used for inter- 
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FIG. 1: Wireframe of a typical considered geometry. Adopt- 
ing the capacitance model we assume that gate, leads and 
nanoparticles are ideal conductors while the other parts of 
the system are modelled as dielectrics. Gate and leads are 
quarter ellipsoids. The substrate is a dielectric cuboid with a 
relative permittivity of 3.9 (silicon dioxide). For the nanopar- 
ticles we assume a spherical geometry. The ligand shell which 
is necessary to bind the nanoparticles to the DNA is modelled 
as a concentric dielectric shell of 0.3nm thickness with a rel- 
ative permittivity of 3. This is shown in the inset of Fig. 
the metallic nanoparticle itself is depicted in grey while the 
space between the particle and the outer sphere is filled with 
the dielectric. Finally the DNA is modelled as a dielectric 
cylinder with a diameter of 2nm and a relative permittivity 
of 3. (The figure was made with "FastCap" 9 .) 

esting electronic applications 17 . Typically the IV charac- 
teristics for an array with Z nanoparticles have the form: 

I(V) oc (V/V T - l) a Q{V-V T ) , V T ozZ b E c (1) 
Here Vr is the threshold voltage, which borders the re- 
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gion of the Coulomb blockade - in a certain bias voltage 
range no current flows due to the Coulomb interaction be- 
tween charged nanoparticles, E c is the charging energy, 
which is a measure for the energy necessary to charge a 
nanoparticle with an additional elementary charge. 

In the studies mentioned above, the electronic trans- 
port is addressed within the so-called orthodox theory. 
Our treatment differs substantially from that theory in 
two aspects: firstly it is not needed to assume that the 
mean occupation of the electronic levels on the nanopar- 
ticles is given by a Fermi distribution. Instead it is in 
general determined both by the tunneling kinetics and 
dissipative effects on the nanoparticles. If the latter dom- 
inate, they drive the mean configuration to a Fermi dis- 
tribution which is the limit used in the orthodox theory. 
But within our model the strength of the dissipation can 
also be weak. We focussed on the latter case in which 
the mean occupation on a nanoparticle is determined by 
the tunneling kinetics alone. Secondly, in the orthodox 
theory the electronic spectrum is assumed to be contin- 
uous. In this paper, however, the nanoparticles are con- 
sidered to be so small that the level spacing can have 
the same order of magnitude as the thermal energy ksT 
or the charging energy E c . So we take into account ex- 
plicitly the discrete nature of the electronic spectrum of 
the nanoparticles. Doing so, the number of many-particle 
states that may take part in transport is so huge that the 
problem completely defies an analytical solution. Instead 
we use a Monte Carlo method to determine the electron 
transport properties. The algorithm used here is differ- 
ent from the one used in previous studie a 10 ! 11 ' 16 : it copes 
with the hugeness of the state space and is partially self- 
adapting to the size of the examined system. As in the 
orthodox theory, tunneling is treated as a perturbation 
while the Coulomb interaction between charged nanopar- 
ticles is taken into account nonperturbatively within a 
capacitance model. In this model the capacitance matrix 
determines the classical effects of Coulomb interaction 
and therefore the transport properties of the system. To 
elucidate the relation between the shape of the device 
and its transport properties the capacitance matrix is 
extracted numerically from realistic array geometries. 

With our model we are able to identify one-electron 
levels and study the interplay of the charging energy and 
the one-electron level spacing. We will further demon- 
strate the impact of dissipation on the IV characteris- 
tics. We will present several effects which are related to 
the geometry and are therefore designable: the IV char- 
acteristics are strongly asymmetric due to asymmetric 
capacitance matrices and the level spacing varying over 
the array. We will show that in the investigated geome- 
try the conductance of the ohmic parts of the IV curves 
may unexpectedly rise with a growing array length. A 
designable NDC effect occurs if finite electronic spectra 
on the nanoparticles are considered. 



I. MODEL 
A. Model system 

In Fig. 0a pictorial representation of the model system 
is shown. The system is modelled by a tight-binding tun- 
neling Hamiltonian for spinless electrons. The free part 
consists of the electronic spectra of the reservoirs and the 
nanoparticles as well as the classical Coulomb interac- 
tion. The reservoirs shall consist of non-interacting elec- 
trons with a continuous, homogenous, infinite spectrum 
and they shall be in thermal equilibrium at all times. 
Their occupation numbers accordingly obey a Fermi dis- 
tribution. The one-electron spectra on the Z nanoparti- 
cles are considered to be discrete and on dot i we con- 
sider explicitly Zi levels with level spacing Ae*. The 
Coulomb interaction is incorporated by the capacitance 
model which is detailed in the next section. We consider 
mutual capacitances between nearest neighbours like e.g. 
C12 as well as between distant conductors like e.g. Cla- 
Potentials can be applied to the reservoirs &r) and 
a gate (&g)- The uncontrollable influence of the DNA 
creates background charges on the dots (indicated by the 
random potentials $> g i). 

The perturbation comprises three types of transitions: 
Transitions within the array We consider phonon as- 
sisted tunneling between the one-electron levels of nearest 
neighboring dots (w TA ). The tunneling matrix element 
t A is assumed to be equal for all tunneling transitions 
within the array. The phonon assisted tunneling stems 
from a linear coupling of a bosonic bath (bosonic bath 2) 
and the electronic degrees of freedom in the arraj*i2i2S. 
Here the bosonic bath shall model phonons in the sub- 
strate. 

Transitions between array and reservoir We include 
tunneling between the continuous reservoirs and the one- 
electron levels of the outermost nanoparticles where no 
phonon assistance is considered here (w TRes ). The tun- 
neling matrix element t Res for these transitions is consid- 
ered to be equal for both reservoirs but generally different 
from t . 

Transitions on a single dot These transitions among 
the one-electron levels on a single dot (w rel ) can be jus- 
tified microscopically by a Frohlich-Hamiltoniani^ (cou- 
pling to bosonic bath 1). 

We assume that the tunneling matrix elements have 
phases which fluctuate randomly due to slight temporary 
changes of the array geometry. Implicit summation over 
these phases prohibits any first order contributions of the 
tunneling matrix elements to the perturbation expansion. 

B. Transition energies 

The transition rates which are given in the next sec- 
tion are determined by the transition energies, i.e. the 
changes in energy of the array that accompany the tran- 
sitions of an electron from one level to another. It is 
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FIG. 2: Model system 



composed of a change in electrostatic and one-electron 
energy. The former is computed within the capacitance 
model: the reservoirs, the gate and the dots are treated 
as macroscopic conductors with potentials $^ and total 
charges Qi. We call the gate and the reservoirs volt- 
age nodes (with potentials J?y and charges Q ) since 
their potential is fixed, the dots are called charge nodes 
(with potentials $ c and charges Q c ) since their charge is 
known and their potentials have to be determined^. The 
potentials and the charges are related by the capacitance 
matrix C 



Q 



c ) =C 

V 



Cb C a 



(2) 



The relation between the mutual capacitances Cy indi- 
cated in Fig.|21and the elements of the capacitance matrix 



Cij is given by Cy = £y J2k=a C jk 
static energy can be written as: 



Cm. The electro- 



E el = ^ (q_ c + & c + qS) t c- 1 (q c + q_ c + q*) 

(3) 

where Q' c := — § v is the polarization charge induced 

by the potentials $ y on the voltage nodes and Q b ^ are 
fixed background charges which account for the electro- 
static influence of the DNA and unknown fabricational 
details. The potentials on the charge nodes are then: 



(4) 



Furthermore we have to take into account the electronic 
spectra of the nanoparticles: the one-electron energy en 



of level I on dot i is En = epi + I ■ Ae^ with the Fermi 
energy epi on dot i. We find the following transition 
energies^: 

Transitions within the array from level / on dot i to 
level on dot i ± 1 

AE ±,i,l'l = y . A£ . ±i _ l . A£j 

-e (4> Cj± i - $c0 + y (C^ 1 - 2Cr i ^ 1 + C~ ± \ l±1 ) 

(5) 

Note that the dependence of the potentials $ c on the 
charges is implied and that we include the constant 
terms EFi m the definition of J> c , i.e. we treat them like 
background charges: 

*c - C- 1 (q_ c + Q; c + <£ 9 ) with <g = Q b ^ - -C c e F 

(6) 

where the vector e F contains the constants SFi- 

Transitions between array and reservoir from/to lead 
a into/out of level I on the neighbouring dot 

AE ±, a ,i = ±l . A£0 ± ( _ e) ^ _ + c -i?_ (7) 

where (3 — 1 for a — L (left lead) and f3 = Z for a = R 
(right lead). (Note that the nanoparticles are numbered 
consecutively from left to right.) 

Transitions on a single dot from level I to level V on 
dot i only change the one-electron energy: 



AE h 1 = {V - I) ■ A £i 



(8) 



It is convenient to express the transition energies in this 
way since we have to calculate the potentials $ c only 
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once at the beginning of the simulation with © and then 
update them after each transition which is computation- 
ally cheap. 

Obviously, the behaviour of the system depends 
strongly both on the values of the mutual capacitances 21 
and the one-electron level spacings. Concerning the for- 
mer we determine the capacitance matrix numerically 
with the help of "FastCap"- using the geometry shown in 
Fig. ^ Only qualitatively we estimate the magnitude of 
the one-electron level spacing of a single, neutral, isolated 
nanoparticle^: 

A£l = . J_ = C c ^ Qj 3Q _ 1Cr 28 el/m 3 

2m e (37r 2 n i ) 3 r i r * 

(9) 

where nt is the electron density, n the nanoparticle radius 
and the numerical value for c was calculated for bulk gold. 
DFT calculations^ and experiments^ suggest that this 
estimation can at least reproduce the correct order of 
magnitude. 



II. METHOD 
A. Transport theory 

We assume that the tunneling rates defined below 
are much smaller than (fcsT)/S so that we can treat 
the tunneling part of the Hamiltonian as a perturbation 
(weak coupling regime) and consider only the lowest non- 
vanishing order of the perturbation expansion (sequen- 
tial tunneling regime). The reservoir and bath degrees 
of freedom of the density matrix are traced out which 
results in the reduced density matrix. We neglect its 
non-diagonal elements, which is justified if the broaden- 
ing of the levels due to tunneling is small compared to 
the level spacing. The diagonal elements P s can then 
be interpreted as the probabilities of the array states |s) 
where the array state s is given by all occupation num- 
bers nu of level I on dot i; \s) = \{nu} u ). The rates of 
electron transfer are calculated in golden rule approxima- 
tion, thereby assuming that transport is incoherent and 
no coherent eigenstates are formed which stretch over 
the whole array, comparable to molecular orbitals. This 
is justified since in reality there are certainly processes 
which destroy the phase coherence - e.g. slight tempo- 
rary changes in the geometry of the array. 

Under the given assumptions the probabilities P s obey 
a master equation in the stationary limit 25 

P S =J2 ( W ss'Ps> - W S , S P S ) = (10) 
s' 

with golden rule rates itVs from array state s to array 
state s'. Each transition rate belongs to exactly one of 
the following sets: 

Transitions within the array from level I on dot i to 



level I' on dot i ± 1 

W TA s w ±,i,l',l = T Ap^ AE ±,i,l>,^ ,T A = ^- \t A \ 2 

(ii) 

where t A is the tunnneling matrix element within the 
array. The function P(E) is the probability per en- 
ergy for the exchange of energy E with the bosonic 
bath and therefore has to be normalized and must ful- 
fill the condition of detailed balance^ f_ dEP(E) = 

1 , P(-E) = e~ 0E P(E). The latter property stems 
from the nature of the bosonic bath: while it is always 
possible to emit energy, there must be excited bosons in 
the bath if energy shall be absorbed. The P(E) function 
given above is an assumption which is preferably simple 
and has the required properties. 

Transitions between array and reservoir from/to lead 
a into/out of level I on the neighbouring dot 

w TRes = w ±,a,t = pa f ± ( ±A £±,«.l) , T a = — \t Res \ 2 d ( 

(12) 

where t Res is the tunnneling matrix element between a 
reservoir and the neighbouring dot, d a is the density of 
states in the reservoir a and f + (E) = (e° E + 1) , f~ = 

Transitions on a single dot from level I to level V on 
dot i 

w rei = w i,l'i = Y rel (Q(/\E l ' l ' l )g(AE l > Vl ) 

+ <d(-AE l ' l ' l ){l + g(-AE i,l>1 ))) 

(13) 

where T rel is the inverse relaxation time, g(E) = 
(e^ E — 1) _1 and 0(E) is the Heaviside step function. 
Since, as in the case of tunneling within the array, 
these transitions are possible due to the coupling to a 
bosonic bath, the rates fulfill detailed balance: w 1,1 1 = 
e -0AE'- 1 w i,u' f or AE 1 ' 1 ' 1 > 0. If the nanoparticles were 
isolated, i.e. if there were no transitions among them, 
it follows from the detailed balance property that in the 
stationary limit the occupation numbers nu obey a Fermi 
distribution. So while the tunneling transitions drive the 
electron distribution out of equilibrium the transitions on 
single dots effectively cool the electrons. 

B. Monte Carlo algorithm 

If we consider Zi levels on dot i there are pos- 
sible array states. So the master equation ifTUI) defies a 
direct numerical solution except for very small systems. 
Therefore we employ a Monte Carlo (MC) method to re- 
trieve the quantities of interest: the current and the shot 
noise. The key idea of the MC method 2 ^ in this context is 
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to discretize time and to get from the transition rate w a > s 
a transition probability tt s i s by multiplying the rate with 
a finite time step At(s) which may depend on the present 
system state s: 7r s / s (At(s)) = w s ' s At(s) The time step 
has to be chosen sufficiently small so that for the total 
probability 7r s (Ai(s)) to leave the state s it holds 



7r.(Ai(s)) = ^7T s , s (Ai(s)) < 1 Vs 



(14) 



The transitions among the array states, which are gov- 
erned by the probabilities 7iy s , constitute a stochastic 
process which can be simulated with the help of ran- 
dom numbers. In contrast to the MC method used in 
statistical physics^ which samples a (grand) canonical 
ensemble, the system examined here is generally out of 
equilibrium. In each valid MC algorithm the probability 
7r s (At(s)) must be properly represented. We write it as 
follows: 



,(At(a)) = uv.At(a) 



1 



w D(s) 



D(s)w At(s) 



(15) 

where D(s) is the number of possible transitions out of a 
state s and wo is an upper bound for the transition rates 
which is also called attempt frequency. If we choose the 
time step At(s) as At(s) = l/(D(s)wo), the transition 
probability reduces to: 



■K s . s (At(s)) 



1 



w D(s) 



(16) 



Note that with the upper choice of the time step the re- 
quirement 114f) is fulfilled since w S ' S < wo and the number 
of nonzero addends is D(s). It is straightforward to imag- 
ine the corresponding algorithm if you regard each factor 
in (I16|l as an independent probability, see Fig. [21 We 
found that even though our algorithm needs two random 
numbers per time step and an attempt transition may 
be rejected, see step 3 in Fig. [31 it is convenient since 
only one rate has to be computed per time step. For 
systems with many possible transitions it is faster than 
algorithms which need the transition rates of all possible 
transitions in each time stepiSiii. Furthermore our al- 
gorithm has the convenient property to be self-adapting 
to the system size: since Ai(s) oc (D(s))~ 1 the number 
of performed steps scales with the amount of possible 
transitions if the runtime of the simulation is fixed. 

Despite the optimized algorithm we can take only a fi- 
nite number of one-electron levels into account, of course. 
In most cases, however, we want to model wide bands 
on the dots, i.e. we do not want to observe any impact 
of the finiteness of the electron spectrum. Practically 
we increase the number of considered levels until the re- 
sults become steady for the maximum bias voltage which 
we want to apply. In order to minimize the computing 
time, i.e. real time, we consider as few levels as possible, 
of course. Therefore we center the spectra around the 
highest occupied levels in the ground state (i.e. the state 
assumed in equilibrium at zero temperature) since for 
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i. Equilibrate until t = T equ i then set t := 

ii. Sample until t = T meas 

FIG. 3: MC algorithm: 1. Augment time by an amount of 
At(s) = 1/ (woD(s)), 2. Choose one of the -D(s) transitions 
with equal probability 1/D(s), in other words choose an at- 
tempt configuration the system might transit to, 3. Accept 
the choice with probability w s i s /wo, 4. If the transition is 
accepted, update the state of the system, especially the total 
number of possible transitions D(s). 



small bias voltages and low temperatures only low-lying 
excitations about the Fermi edge appear. To determine 
the ground state the total energy of the system which 
includes the electrostatic energy © and the one-electron 
energies has to be minimized. Due to the discreteness of 
the charges Q c this is a non-trivial minimization problem 
which we do not address here. 

The initial state of a simulation run is always the 
ground state, i.e. all spectra are half-filled. Each run 
starts with an equilibration period T equ i in which we let 
the system evolve, by iterating the algorithm (Fig. |3J, 
without sampling the assumed states. So the system can 
reach its stationary state. In the subsequent measure- 
ment period T meas the quantities of interest are retrieved. 
The current can simply be obtained by counting the 
electrons that are transferred e.g. between the left lead 
and the first nanoparticle and dividing by the length of 
the measurement period T meas : I L = Qi(T mea s) /T meas 
where Ql (T meas ) is the charge that is transferred during 
the measurement time T meas . Note that in the station- 
ary state the current through all tunneling barriers is 
the same. For a fixed set of parameters the simulation 
is repeated with different seeds for the random number 
generator in order to get statistically independent runs. 
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With this ensemble we can determine the statistical stan- 
dard error of the mean current. The (zero- frequency) 
shot noise Si L (0) and the Fano factor F can also be 
estimated^: 



the current as 



Sl L (0) = 7^ {(QdTmeasf) - (Q L {T meas )) 2 ) 

J- meas > ' 



F 



2e (II) 



(17) 



where (...) denotes the ensemble average. As for the 
current, in the stationary limit the shot noise is the same 
for all tunneling barriers. Note that all given quantities 
are estimators which become exact in the limit T meas — > 
oo. 

The validity of our method was checked by compar- 
ing our results with the solution of the master equation 
for a small system. With the same benchmark sensi- 
ble values for the simulation parameters (equilibration 
and measurement time, size of the ensemble) were ob- 
tained. Due to the self-adaptive property of the algo- 
rithm these parameters are also suitable for bigger sys- 
tems. For each geometry we increased the number of 
considered one-electron levels until the IV curves became 
steady. Furthermore we found that the computing time 
scales linearly with the number of nanoparticles Z and 
quadratically with the maximum number of considered 
levels on a dot (maxl-^},). It is exponentially smaller 
than the computing time needed for the direct solution 
of the master equation which scales as (2^i Zi ) z . How- 
ever, note that the MC method is not equivalent to a 
solution of the master equation: though the current and 
shot noise may be retrieved from a MC simulation, it is 
practically impossible to determine all probabilities P s 
correctly. 



C. Charge states 

To interpret the simulated results we draw a sample 
out of a single simulation run and look at the proba- 
bilities of charge configurations and mean rates among 
them. To get the probability Pc of a certain charge con- 
figuration C = {Qci}i, we sum the MC times during 
which this configuration is assumed and divide by the 
total MC time. The mean rate w^Iq through a tun- 
neling junction i from state C to state C to the right 
or left respectively is determined by counting the transi- 
tions between state C and state C by tunneling in the 
given direction through junction i and dividing the sum 
by the MC time that is spend in the state C. The cur- 
rent through a tunneling junction i can then be written 
as Ii = J2c c So- a " w c'C^'c with er = ±. In the sta- 
tionary state the currents through the junctions are all 
equal to one another Ii = L Vi = j, so we can write 



I = 



1 



1 ^ 



h.= 



1 ^ 

C,C> 



Pc'Wc'C 



c,c 



. aw 



C'C 



and I, 



C'C 



(18) 

z^Pcwcc- The 



with w C 'c = J2i 

partial current Ice from charge state C to charge state 
C has a positive sign if it flows from left to right and 
the opposite sign for the opposite direction. If there are 
partial currents which flow between the same states - 
then they have necessarily opposite directions, we keep 
only the difference of them so that there is only one net 
partial current between two charge states. 



III. RESULTS 

The following conventions hold for all shown results. 
The general geometrical setup is the one already shown 
in Fig. ^ All energies are normalized to the maximum 
level spacing that occurs in the array (max {Aerj}). In- 
stead of voltages (potentials) we use potential energies 
eV and charges are given in units of e. We open the bias 
voltage window symmetrically (i.e. $l = —Qr) because 
we do not want to introduce artificially an additional 
asymmetry. We define V* = (e&n)/(Ae max ) so that the 
current is positive (i.e. flows from left to right) if V* is 
positive. We divide all calculated golden rule rates by 
T A ■ max(P(A£')), so that the maximum rate within the 
array is equal to 1. We set the maximum rate between 
array and reservoir equal to 0.1 since we assume that the 
tunnel coupling within the array is stronger. Current and 
charge are expressed in units of the elementary charge e. 
We normalize the current to the maximum rate in the 
bulk of the array 7* = (I/e)/(T A ■ xaax(P(E))). The 
current is defined to be positive if it flows from left to 
right. The parameter V of the function P(E) is set to 2, 
which is small enough to see individual one-electron lev- 
els and big enough to give a sufficiently high current. No 
error bars appear in the following results since the rela- 
tive statistical error is always only about 0.1%, so error 
bars would not be visible. 



A. Generalizations of results for single quantum 
dots 

1. Interplay of one- electron level spacing and charging 
energy 

In a 2 nanoparticle array we investigate the case Aei < 
|^-(C _1 )ii|, i.e. the level spacing is smaller than the typi- 
cal charging energy. In agreement with results for a single 
do& we find that the level spacing imposes a fine struc- 
ture on the Coulomb staircase which is related to the 
charging energy. In Fig. ^ the corresponding IV char- 
acteristics are shown. Looking at the partial currents 
Ic'c defined above we find that for the first 4 steps there 
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FIG. 4: IV characteristics of an array of 2 nanoparticle with 
diameters of lnm and \.2nm respectively. The level spacing 
of the first equals about 28007^ so that a voltage V* of 1 
equals about 0.24V . In the inset V* ranges from —10 to 
10. The region in the dashed rectangle is shown in the big 
plot: V* ranges from —3 to +3 and the IV curve is shown 
for T = and T = 0.05 . The gate voltage was tuned so 
that the induced charges on the dots are positive and slightly 
smaller than an integer value in order to have a small Coulomb 
blockade region. 

is only one relevant transport path in the charge con- 
figuration space: (0,0) -> (-1,0) -» (-1,+1) -> (0,0) 
where the charges are given as differences to the ground 
state charges. So the first 4 steps must be due to the 
level spacing. On the 5. step a second path is relevant: 
(0,0) -> (0,+l) -» (-1.+1) -» (0,0). For higher tem- 
peratures the fine structure - i.e. the features due to the 
one-electron levels - is smeared out, see Fig.QJ while the 
typical Coulomb staircase remains: in the middle of each 
plateau a new transport path becomes available. 
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FIG. 5: Array with 3 nanoparticles, comparison of IV curves 
for different relaxation rates at a fixed gate voltage. The 
curves are obtained for values of the relaxation rate prefactor 
F re! between (highest curve) and 20 (lowest curve) and since 
the change of the curves' shape is systematic, they are not 
distinguished. 




FIG. 6: Array geometry 



That is the reason why the steps in the IV curve become 
distinct. Electrons above the Fermi energy relax to lower 
lying levels and the corresponding transition energy is 
dissipated in the bosonic bath. Such electrons can per- 
form fewer transitions as without dissipation since less 
energy is available. That is the reason why the current 
generally decreases with an increasing relaxation rate. 



2. Influence of dissipation 

We consider a 3 nanoparticle array at T = and fixed 
gate voltage and study the impact of a finite relaxation 
rate, see Fig. Generalising results for a single dot 30 
we find that without relaxation the electrons overheat 
and consequently the structures in the IV characteris- 
tics are smoothed. Strong relaxation, on the other hand, 
which effectively cools the electrons, sharpens the steps 
and decreases the absolute value of the current. The 
IV curves for intermediate relaxation rates lie between 
the curves belonging to the extreme cases. These ten- 
dencies can also be observed for other gate voltages and 
other array lengths. They can be understood by noting 
that a high relaxation rate keeps the mean occupation 
of the one-electron levels on a dot close to the equilib- 
rium i.e. Fermi distribution. For T = this leads to the 
formation of a defined Fermi edge on the nanoparticles. 



B. Results uniquely related to the array geometry, 
designable effects 

We examine arrays of nanoparticles with uniformly 
growing diameters, see Fig.Elfor the case of 4 nanoparti- 
cles. We choose this special geometry for two reasons: on 
the one hand the array can be enlarged in a defined way. 
We start with two nanoparticles and let the number of 
nanoparticles and therefore the array length grow so that 
we discover certain features which evolve systematically 
with increased length. On the other hand it is interesting 
to combine small and big nanoparticles since they differ 
strongly both in level spacing and charging energy. In 
the case of 5 nanoparticles, which is the longest array 
that is studied here, the diameters of the nanoparticles 
range from lnm to 1.8nm in steps of 0.2nm. The level 
spacing of the smallest particle is estimated to be about 
Q.2AeV according to eq. iJjJJ so that a voltage V* of 1 
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I* 




FIG. 7: Array of 4 nanoparticles, comparison of IV curves for 
different gate voltages. Upper Linear fit of the ohmic part of 
the curves (dashed grey line) has the slope G which is called 
overall conductance (denoted by the slope of the small trian- 
gle). The 1/*-intercept of the fitted line is the offset voltage 
V ff (denoted by the double-headed arrow). Lower The cur- 
rent is multiplied by a factor and the curves are artificially off- 
set from each other proportional to the applied gate voltage. 
The flat region in the middle always corresponds to Coulomb 
blockade, i.e. 7* = 0. The dashed line borders the Coulomb 
blockade region as a guide for the eye. 



in the following curves equals then 0.241/. This energy 
equals a temperature of about 2800-ftf so that we find 
level spacing related features in the current or shot noise 
at finite temperatures. Since the size of the nanoparticles 
increases from left to right, the level spacing decreases in 
the same direction. 



1. Strong asymmetry of IV characteristics 

In Fig. [7| the IV characteristics of an array of 4 
nanoparticles exemplify the typical IV characteristics of 
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TABLE I: Offset voltage V g and overall conductance G with 
respect to the number of nanoparticles Z 



the array geometry. The various curves correspond to 
different gate voltages. We observe that for all array 
lengths two regions with different behaviours can be dis- 
tinguished. For smaller bias voltages the curves exhibit 
a striking asymmetry with respect to reversal of the bias 
voltage. For example the threshold voltage is different for 
positive or negative bias voltages (V^r). This asymmetry 
is more pronounced for longer arrays and for curves with 
higher threshold voltages, especially if the V^r exceed 
the offset voltage V g. The position, width and height 
of the steps in this bias voltage region and V^r depend 
strongly on the gate voltage. For larger bias voltages, 
however, the IV curves shown in Fig. approximately 
coincide and become symmetric (with respect to reversal 
of the bias voltage). For small bias voltages the number of 
many-particle states or the number of paths through the 
charge state space that take part in transport is smaller. 
So it is important which states or paths actually partici- 
pate. This is in turn influenced by the gate voltage since 
it shifts the energetic position of the charge states and 
determines therefore which paths are available. For big- 
ger bias voltages many states or paths participate so it 
should be less important whether a certain state takes 
part or not: what we observe is their "average" contribu- 
tion. That is why the curves for different gate voltages 
coincide for high bias voltages. It is also the reason for 
the disappearance of the asymmetry in the same region. 
A detailed discussion of the asymmetry can be found in 
the appendix B 



2. Overall conductance and offset voltage 

As already mentioned above, the IV curves shown in 
Fig. \7\ approximately coincide and become symmetric for 
larger bias voltages. In that region the current increases 
linearly, apart from very small steps. We fit this linear 
segment with a straight line and define its slope as the 
overall conductance G. The offset voltage V g, i.e. the 
V"*-intercept of the fitted line, can be thought of as a 
kind of "mean" threshold voltage. The actual threshold 
voltage of a single curve obviously depends on the gate 
voltage, as already mentioned. G and V g are approxi- 
mately the same for positive and negative bias voltages. 
Our results for G and V g with respect to the number 
of nanoparticles Z are compiled in the two left columns 
of Table [I] for the case without relaxation (w = 0) . 
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We find that G and V g both increase with an increasing 
array length. 

The increase of V g with an increasing array length 
has also been observed elsewher o 15 ' 16 so it seems to be 
quite generic for arrays. One factor which contributes 
to this tendency is the decrease of the capacitances with 
increasing distance between the conductors. Especially 
the capacitances that couple the leads with the nanopar- 
ticles C a i, see eq. (|2J|, are important. The cpacitances ex- 
tracted from the geometry with the help of "FastCap" 9 
do not decrease linearly, like e.g. for a simple parallel 
plate capacitor, but approximately exponential. This is 
reasonable since the nanoparticles partially screen the 
electric field. For two neighbouring nanoparticles that 
are in the middle of the array the difference between the 
coupling capacitances is small. So we have to apply a 
high bias voltage to create a potential difference between 
these two particles which permits a tunneling transition, 
see eq. <|11[) . Therefore the threshold voltages and V g 
increase with increasing array length. Concerning G the 
change of the number of transport paths with increas- 
ing bias voltage plays an important role. This change 
is in turn determined by the change of the number of 
possible transitions. A certain transition between two 
neighbouring dots becomes possible when the potential 
gradient between the dots becomes high enough. If we 
assume that this gradient is roughly the bias voltage di- 
vided by the number of tunnel junctions then we should 
expect that G decreases with increasing array length. On 
the other hand we "grow" the array by adding bigger 
nanoparticles with a higher density of one-electron states. 
Given that the number of extra electrons in the array is 
small a higher number of one-electron levels within the 
bias voltage window results in a higher current. This 
effect might overcompensate the reduced potential gra- 
dient which would result in the observed behaviour of G. 
To check this assertion we have artificially raised the level 
spacing on the last dot in the 3 particle array so that it is 
equal to the level spacing on the first dot: in this case we 
find that the overall conductance G is only 0.1 compared 
to 0.14 with the original level spacing. So G can indeed 
be lowered by raising the level spacing on the last dot. 
This supports our assertion. Obviously the offset voltage 
and the overall conductance can be tailored by the choice 
of the nanoparticle sizes. 

For the case of high relaxation rates G and V g are 
compiled in the two right columns of TableQ] Comparing 
with the case without relaxation we observe that the off- 
set voltage is generally bigger if there is relaxation while 
the overall conductance is approximately the same. As 
argued in section ITlI A 21 the current generally decreases 
with an increasing relaxation rate. This homogeneous 
downward shift of the IV curve correspondingly increases 
the offset voltage while the overall conductance remains 
unaltered. 
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FIG. 8: 4 nanoparticle array, striking asymmetry in the IV 
characteristics (upper) which is accompanied by a giant Fano 
factor (lower). The temperature of T = 0.1, for which also 
the Fano factor in the lower figure is calculated, equals about 
2807V 



3. Giant step, giant Fano factor 

As already stated above we observe that in the IV char- 
acteristics for several gate voltages the first current step 
to the left or to the right of the Coulomb blockade re- 
gion is strikingly high. One example is shown in Fig. [S] 
Furthermore the higher step is accompanied by a peaked 
giant Fano factor. To determine the origin of the big step 
in the current and the giant peak in the Fano factor we 
record the dynamics of the simulation, i.e. in Fig. the 
transferred charge is plotted against the MC time for a 
single simulation run. This approach has already been 
used in a MC study of the electron transport properties 
in a different model^i in which also a giant Fano fac- 
tor is found. In Fig. we have recorded the dynamics 
for the voltages marked in Fig. [5] For the top of the 
peak marked with 2. we observe comparatively long pe- 
riods without charge transfer and intermediate tunneling 
avalanches. Looking at the definition of the Fano factor 
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5000 



FIG. 9: Transferred charge vs. time showing tunneling 
avalanches, bias voltages correspond to those marked in 
Fig. Hi. V* = -2.5, 2. -> V* = 4.5, 3. ->T = 5.0, 4. 
-> V* = 5.5 



Current carrying states 




Quasi blocking state 



FIG. 10: Quasi-blocking state, the width of the arrows cor- 
responds to the size of the transition rates. A quasi blocking 
state is characterized by numerous, big ingoing rates and few, 
small outgoing ones. In our case, the rates among the states 
that are involved in the charge transport are also much bigger 
than the rate leading out of the quasi blocking state. 



F, eq. Q17|). there are two ways to conclude that the ob- 
served dynamics result in a super poissonian noise (i.e. 
F > 1). On the one hand, one can regard the avalanches 
on a much bigger time scale as single, statistically inde- 
pendent transfer processes in which an effective charge 
bigger than e is transported. Then one uses the Poisson 
value for the shot noise 2q(I) with an effective charge 
q bigger than e which results in a Fano factor F > 1. 
On the other hand, one can regard the avalanches as 
a bunching of the tunneling processes and this positive 
correlation leads per definition to F > 1. So both inter- 
pretations deduce the super poissonian shot noise from 
the observed dynamics. For bigger bias voltages, 3. and 
4., the length of the periods without charge transfer is 
reduced so the Fano factor is reduced correspondingly. 
The big step height can be understood analogously: dur- 
ing the periods without charge transfer a quasi-blocking 
state is assumed , i.e. the sum of rates that lead out of 
this state is much smaller than the sum of the rates that 
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V* 



10 20 30 40 y* 

FIG. 11: 2 nanoparticle array, NDC effect is stronger for 
higher tunneling rates between array and reservoirs (w TRes ) 
(upper) whereas the NDC effect is attenuated by a broad 
P(E) function i.e. for a large F (lower) 

lead into the same state, see Fig. ^3 So if the state is 
visited the system rests there for a comparatively long 
time. If no rates led out of that state it would be a real 
blocking state and the current would be zero since the 
dynamics of the system would ultimately end up in this 
state. The probability P s of a blocking state is therefore 
1. In the shown case we find that a neutral many-electron 
state (0,0,0,0) serves as a quasi-blocking state (respec- 
tively blocking state in the region of Coulomb blockade). 
If the system leaves it, current can flow. This current is 
however carried by other states with much higher rates 
among them. Therefore the high step in the IV charac- 
teristics. 



4- Designable NDC effect 

So far we have modelled infinite spectra, see section 
III Bl Now we intentionally consider only a few levels on 
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(P(E) 




FIG. 12: 2 nanoparticle array, the energies {I ■ Aet + e$ci} ; 4 
are sketched. The P(E) function determines the transition 
rates between the sketched levels. Upper: low bias voltage 
and small F. Middle: high bias voltage and small V, the 
high bias voltage results in a large charge gradient between 
the nanoparticles which causes the large offset between the 
spectra. Lower: high bias voltage and big F, a broad P(E) 
function attenuates the NDC effect. 



the higher the rate of tunneling between array and reser- 
voirs (w TRes ). So this is an example of an effect which 
strongly depends on the geometry and which is therefore 
designable. The origin of the NDC can be understood 
by looking at the transition energies, see section llBl In 
Fig. El w e sketch the energies {I ■ Aei + e$ci}/ i f° r the 
two nanoparticles. For T = a transition between the 
nanoparticles can only happen if the difference between 
those energies is bigger than (C^ 1 - 2(7^ + C~± u±1 ), 
see eq. (JSJ . For finite temperatures this restriction is soft- 
ened by the coupling to the bosonic bath which results 
in the P(E) function. The offset between the sketched 
spectra is determined by the difference between the po- 
tentials e$ci which depend in turn on the charges on 
the nanoparticles. If the tunneling rates between array 
and reservoirs are sufficiently large a high charge gradi- 
ent accumulates and therefore a large offset between the 
spectra occurs. Since the spectra are finite, the transi- 
tion energies of possible transitions are big if the offset 
between the spectra is large. The corresponding value of 
the P(E) function is small. Consequently the transition 
rates among the nanoparticles become smaller with in- 
creased voltage which results in the NDC. If the P(E) 
function is broad, this effect is softened, see eq. JTJl. 
The higher the factor T the broader the P(E) function. 
Therefore the NDC effect is less distinct for higher values 
ofT. 



IV. CONCLUSION 

We have shown that an improved MC algorithm can 
cope with the complexity of the electronic transport 
through nanoparticle arrays with discrete electronic spec- 
tra. Though we have considered linear arrays only we 
found a huge variety of transport functionalities. Of 
course, other geometries would also be worth studying, 
e.g. ring-shaped devices^ can be used as charge storage 
elements. Since arrays can be made self-assembling, are 
robust against fabricational imperfections and show, as 
we found out, transport features that can be designed 
via the geometry, we consider them to be ideal building 
blocks for electronic architecture on the nanoscale. 



each nanoparticle. This is interesting since the electronic 
structure of small nanoparticles differs in general strongly 
from bulk material^. We assume on each nanoparticle 
a finite "band" of discrete, equally spaced energy levels. 
The level spacing is estimated as given above, see eq. (0 ■ 
Under these assumptions we find for an array of two 
nanoparticles a NDC which depends strongly on the ratio 
between transition rates within the array (w TA ) and the 
rates between array and reservoirs (w TRes ), see Fig. 1111 
This ratio can be tuned by varying the distance between 
array and contacting reservoirs: the smaller the distance 
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APPENDIX: ASYMMETRY CONCERNING 
REVERSAL OF BIAS VOLTAGE 



The occurence of IV curves which are asymmetric con- 
cerning the reversal of the bias voltage would be surpris- 
ing if the transport was calculated within the Landauer- 
Biittiker formalism 33 - 34 . Within this approach the IV 
curves must be symmetric due to time reversal symme- 
try: for each wave which travels through the system from 
left to right there is a time reversed one with the opposite 
direction of propagation. But the Landaucr-Biittiker for- 
malism can only be applied for coherent transport with- 
out interaction. Both suppositions are violated in our 
case: we assume that the phase information is lost at 
every tunneling event and we include Coulomb interac- 
tion, see sections II Al and 111 Al Therefore, in our case, 
asymmetry is present in general unless special symme- 
tries impose symmetric curves. We will discuss now two 
symmetries that are actually relevant in our case. By ex- 
plicitly showing how these symmetries are broken by the 
Coulomb interaction and the distribution of the density 
of states (inverse level spacing) the occurence of asym- 
metric IV curves is rendered plausible. We adopt the 
framework of the orthodox theory for these considera- 
tions since the basic mechanisms can be understood by 
looking at charge states alone. 



1. Particle- hole symmetry 



This symmetry means that two paths through the 
charge state space with opposite charges are equivalent. 
Equivalent here means that they appear at the same ab- 
solute value but opposite sign of the bias voltage and 
that the rates for corresponding transitions are equal. 
E.g. for a 2 nanoparticle array this means that if the path 
(-1, -2) -> (-2, -2) -> (-1, -3) ->■ (-1, -2) appears at 
V = x then the path (1,2) -> (2,2) (1,3) (1,2) 
appears at V = —x and the rates for corresponding 
transitions are equal. This symmetry is present if there 
is no coupling to the gate and if there are no back- 
ground charges. To conclude this we have to look at 
eq. (0 for the potentials J> c on the nanoparticles: $ c = 

C^ 1 (Q + Q' + Q ) . If there is no coupling to the 



gate the polarization charges on the nanoparticles are re- 
versed if the bias voltage is reversed: V — > — V => Q' c — > 
—Q' c - If we reverse the charges Q c at the same time 

and there are no background charges Q^, the potentials 

% $ c -> -to- 



<& c are reversed: V — » — V, ^ ( , 



Since the potentials $ c determine the transition rates 
two charge states with opposite charges are equivalent in 
the sense explained above. Obviously hnite background 
charges or coupling to the gate break this symmetry. 



2. Inversion symmetry 

This symmetry means that two paths through the 
charge state space which are mirror images of each other 
are equivalent where the mirror plane shall be situated 
in the middle of the array and equivalent is meant in 
the sense explained above. E.g. for a 2 nanoparticle ar- 
ray this means that if the path (—1, —2) — > (—2, —2) — + 
(—1,-3) — > (—1,-2) appears at V = x then the path 
(-2, -1) ->■ (-2, -2) -> (-3, -1) -> (-2, -1) appears at 
V = —x and the rates for corresponding transitions are 
equal. This symmetry is obviously present if the whole 
setup is symmetric with respect to a mirror plane in the 
middle of the array. Both an asymmetric capacitance 
matrix and an asymmetric distribution of the density of 
states on the nanoparticles break this symmetry. 



a. Asymmetric capacitance matrix 

We look at an array of 2 nanoparticles with different 
sizes, the left nanoparticle (A) shall be smaller than the 
other (B). The capacitative coupling (between nodes) 
shall be negligible so that the capacitance matrix is di- 
agonal. Then each nanoparticle can be characterized 
by a single capacitance, Ca and Cb respectively. Since 
nanoparticle A is smaller, it holds Ca <C Cb- The oppo- 
site holds for the charging energies E Ca E Cb so more 
energy is needed to charge nanoparticle A. Now it is clear 
that the paths (0,0) -> (-1,0) -> (0,-1) -> (0,0) and 
(0,0) -> (0,-1) -> (-1,0) -> (0,0) are not equivalent. 
The asymmetric capacitance matrix breaks the inversion 
symmetry. 



b. Asymmetric distribution of the density of states 

The capacitance matrix of the 2 nanoparticle array 
shall now be symmetric and still diagonal, so the charg- 
ing energy is the same for both nanoparticles. The den- 
sity of states, however, shall be Da on the left and Db 
on the right nanoparticle. (This corresponds to differ- 
ent level spacings in our model.) The density of states 
in both reservoirs shall be D Res . The energy change for 
the three tunneling transitions that appear in the path 
(0,0) -> (-1,0) -► (0,-1) -> (0,0) shall be denoted by 
AE\,AE2 and AE$. (The index denotes the number of 
the transition e.g. AEi is the energy difference occuring 
at the transition (0,0) — > (—1,0).) The corresponding 
transition rates are Wi, w 2 and W3. For zero tempera- 
ture the orthodox theory rates are: 



W2 
W3 



2tt 

Y 

2tt 
2tt 



Res I 



,A\ Z 



(-AE 1 )e(-AE 1 )D Res D A (A.l) 
(-AE 2 )0(-AE 2 )D a Db (A.2) 



t Hes {-AE 3 )Q(-AE 3 )D B D Res (A.3) 
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If we reverse the bias voltage and look at the path 
(0,0) -> (0,-1) -► (-1,0) -> (0,0), we get different 
energy differences AE'^AE^ and AE' 3 and rates w[, w' 2 
and w' 3 . (Note that here e.g. AE[ is the energy difference 
occuring at the transition (0, 0) — -> (0, —1).) Since the ca- 
pacitance matrix is symmetric, it holds AE[ = AEi V«. 
So we get: 

w > = ^\t R ^\ 2 {-AE 1 )@{-AE 1 )D Res D B (A.4) 
< = ^|^| 2 (-Ai? 2 )e(-Ai? 2 )^^ B (A.5) 
w , a = ?l\ t Res\ 2 (-AE 3 )Q(-AE 3 )D A D Res (A.6) 



Obviously the rates {wi} i are different from the rates 
{w-},; due to the different density of states {D i } i so the 
considered paths are not equivalent. An asymmetric dis- 
tribution of the densities of states breaks the inversion 
symmetry. 
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